Libraries:

# libraries -------------------
library(dplyr)
library(tidyverse)
library(ggplot2)
library(GGally)
library(glmmTMB)
library(TMB)
library(emmeans)
library(DHARMa)
library(lmtest)
library(ggeffects)
library(performance)

Set working directory

knitr::opts_knit$set(root.dir = '/Users/user/Desktop/Data/Regen_RProj/')

Functions:

# functions -------------------

Source large seedling data:

source("Scripts/LS_Import.R")

Pivot wider to create dataframe where each row is for one plot and has total details for each species group

ls_merge2 <- ls_merge %>% 
  select(Plot_No, Region, Treat_Type, Site, Species_Groups, Total) #this is dropping browse and stump sprout data

ls_merge2 <- ls_merge2 %>% 
  pivot_wider(names_from = Species_Groups, values_from = Total)

Import time since data and add it to the large seedling dataset

time_since <- read_csv("CleanData/Treat_Year_Data.csv")

ls_merge3 <- merge(ls_merge2, time_since, by = 'Site')
#log transform time from treatment data
ls_merge3$l.TFT <- log(ls_merge3$Time_from_Treat)

Run the ‘Add_BA’ script and merge with dataset:

source("Scripts/Add_BA.R")

# merge with ls dataset -------------------
ls_merge4 <- merge(ls_merge3, prism_BA, by = "Plot_No")

Run ‘Ground_Data.R’ script and add it to large seedling dataset:

source("Scripts/Ground_Data.R")

# merge with ls dataset -------------------
ls_merge5 <- merge(ls_merge4, ground3, by = "Plot_No")

Source and add veg data

source("Scripts/Veg_Data.R")

# merge with ls dataset 
ls.all <- merge(ls_merge5, veg3, by = "Plot_No")

rm(ls_merge5,
   ls_merge2,
   ls_merge3,
   ls_merge4)

The large seedling count data is taken in 10m2 plots; basal area is measured in hectares; veg and soil data is taken in 1m2 plots.

Large seedling data will be converted into 1m2 plots in order to compare across and reduce the amount of scales of data collection to two: 1m2 plots and per hectare observations.

Create log transformed categories for newly added variables, then select for just the desired variables:

ls.all$PIRI.1m <- ls.all$PIRI/10
ls.all$SO.1m <- ls.all$Shrub_Oak/10
ls.all$Other.1m <- ls.all$Other/10

ls.all$l.PIRI1 <- log(ls.all$PIRI.1m + 1)
ls.all$l.SO1 <- log(ls.all$SO.1m + 1)
ls.all$l.other1 <- log(ls.all$Other.1m + 1)

ls.all2 <- ls.all %>% 
  select(Treat_Type, Region, Site, Plot_No, PIRI, PIRI.1m, l.PIRI1, Shrub_Oak, SO.1m, l.SO1, Other, Other.1m, l.other1, Time_from_Treat, l.TFT, BA_HA, l.BA_HA, PIRI.BA_HA, l.BA_piri, Mineral_Soil, l.Mineral, Litter_Duff, avgLD, avgLD_l, Veg_Total, l.Veg_Total) %>% 
  arrange(Treat_Type)

Select just for numerical vs log and then look at paired plots:

#not transformed
ls.num <- ls.all2 %>% 
  select(PIRI, Shrub_Oak, Other, Time_from_Treat, BA_HA, PIRI.BA_HA, Mineral_Soil, avgLD, Veg_Total, Treat_Type)

ggpairs(ls.num)
ggpairs(ls.num, aes(color = Treat_Type))


#log transformed
ls.numl <- ls.all2 %>% 
  select(l.PIRI, l.SO, l.other, l.TFT, l.BA_HA, l.BA_piri, l.Mineral, avgLD_l, l.Veg_Total, Treat_Type)

ggpairs(ls.numl)
ggpairs(ls.numl, aes(color = Treat_Type))

rm(ls.num,
   ls.numl)

Can see the correlation coefficients for linear (Pearsons) relationships. None of them appear very strong, except for ones that are analogs (avg LD vs mineral soil exposure; ba/ha vs piri ba/ha)

Log transformed average litter depth and basal area per hectare have a weak relationship (corr 0.33), which does make sense.

The above script a data naming conventions are the same as in the LS_GLMM script

Full dataset is called ls.all2

Starting first with model analysis without treatment type:

Basal area doesn’t seem important, so I’m going to go back to the 10m2 and 1m2 scales and lose per hectare observations. I’ll need to rework the dataset

Revised data set with LS observations at 10m2 scale; this means no non-interger values for the poisson distro

ls.all3 <- ls.all

ls.all3$l.PIRI <- log(ls.all3$PIRI + 1)
ls.all3$l.SO <- log(ls.all3$Shrub_Oak + 1)
ls.all3$l.other <- log(ls.all3$Other + 1)

ls.all3 <- ls.all3 %>% 
  select(Treat_Type, Region, Site, Plot_No, PIRI, l.PIRI, Shrub_Oak, l.SO, Other, l.other, Time_from_Treat, l.TFT, BA_HA, l.BA_HA, PIRI.BA_HA, l.BA_piri, Mineral_Soil, l.Mineral, Litter_Duff, avgLD, avgLD_l, Veg_Total, l.Veg_Total) %>% 
  arrange(Treat_Type)

testing zero inflation

tapply(ls.all3$PIRI, ls.all3$Region, summary)
## $ALB
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.00000 0.00000 0.01639 0.00000 1.00000 
## 
## $LI
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.0137  0.0000  1.0000 
## 
## $MA
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     0.0     0.0     0.0     0.5     0.0    14.0 
## 
## $ME
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0       0       0       0       0       0 
## 
## $NH
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   0.000   0.266   0.000   9.000
# no PIRI in ME
tapply(ls.all3$PIRI, ls.all3$Treat_Type, summary) 
## $Control
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.00000 0.00000 0.04274 0.00000 3.00000 
## 
## $FallRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.00000 0.00000 0.09804 0.00000 4.00000 
## 
## $Harvest
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   0.000   1.304   1.250  14.000 
## 
## $MowRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.1103  0.0000  9.0000 
## 
## $SpringRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.00000 0.00000 0.02299 0.00000 1.00000

I started modeling with poisson distro and excluding treatment type. Models without treatment type failed in model fit. Going to run variable elimination again, starting with models with treatment type

ls.m13 <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + l.BA_HA + l.BA_piri + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = nbinom1)
AIC(ls.m13) #288.9
## [1] 288.9363
ls.m14 <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + l.BA_HA + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = nbinom1)
AIC(ls.m14) #287.8
## [1] 287.8049
ls.m15 <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = nbinom1)
AIC(ls.m15) #285.9
## [1] 285.9361
ls.m16 <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = nbinom1)
AIC(ls.m16) #284.4
## [1] 284.4039
ls.m17 <- glmmTMB(PIRI ~ Treat_Type + l.SO + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = nbinom1)
AIC(ls.m17) #283.9
## [1] 283.9195
lrtest(ls.m13, ls.m14, ls.m15, ls.m16, ls.m17) #non significant
## Likelihood ratio test
## 
## Model 1: PIRI ~ Treat_Type + l.SO + l.other + l.BA_HA + l.BA_piri + l.Mineral + 
##     avgLD_l + l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + l.SO + l.other + l.BA_HA + l.Mineral + avgLD_l + 
##     l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## Model 3: PIRI ~ Treat_Type + l.SO + l.other + l.Mineral + avgLD_l + l.Veg_Total + 
##     offset(l.TFT) + (1 | Site/Plot_No)
## Model 4: PIRI ~ Treat_Type + l.SO + l.other + avgLD_l + l.Veg_Total + 
##     offset(l.TFT) + (1 | Site/Plot_No)
## Model 5: PIRI ~ Treat_Type + l.SO + avgLD_l + l.Veg_Total + offset(l.TFT) + 
##     (1 | Site/Plot_No)
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1  15 -129.47                     
## 2  14 -129.90 -1 0.8686     0.3513
## 3  13 -129.97 -1 0.1312     0.7172
## 4  12 -130.20 -1 0.4678     0.4940
## 5  11 -130.96 -1 1.5156     0.2183
ls.m18 <- glmmTMB(PIRI ~ Treat_Type + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = nbinom1)
AIC(ls.m18) #284.4
## [1] 284.3545
lrtest(ls.m17, ls.m18) #p = 0.12
## Likelihood ratio test
## 
## Model 1: PIRI ~ Treat_Type + l.SO + avgLD_l + l.Veg_Total + offset(l.TFT) + 
##     (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + avgLD_l + l.Veg_Total + offset(l.TFT) + (1 | 
##     Site/Plot_No)
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1  11 -130.96                     
## 2  10 -132.18 -1 2.4349     0.1187
ls.m19 <- glmmTMB(PIRI ~ Treat_Type + avgLD_l + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = nbinom1)
summary(ls.m19) #282.4
##  Family: nbinom1  ( log )
## Formula:          
## PIRI ~ Treat_Type + avgLD_l + offset(l.TFT) + (1 | Site/Plot_No)
## Data: ls.all3
## 
##      AIC      BIC   logLik deviance df.resid 
##    282.4    320.3   -132.2    264.4      489 
## 
## Random effects:
## 
## Conditional model:
##  Groups       Name        Variance  Std.Dev. 
##  Plot_No:Site (Intercept) 7.864e-09 8.868e-05
##  Site         (Intercept) 8.389e+00 2.896e+00
## Number of obs: 498, groups:  Plot_No:Site, 498; Site, 47
## 
## Dispersion parameter for nbinom1 family (): 1.84 
## 
## Conditional model:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -8.0011     1.8849  -4.245 2.19e-05 ***
## Treat_TypeFallRx     2.8946     1.8244   1.587 0.112610    
## Treat_TypeHarvest    5.8075     1.9747   2.941 0.003271 ** 
## Treat_TypeMowRx      2.2997     1.7665   1.302 0.192968    
## Treat_TypeSpringRx   2.9793     1.9726   1.510 0.130956    
## avgLD_l             -0.9013     0.2489  -3.622 0.000293 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
lrtest(ls.m18, ls.m19) #0.85
## Likelihood ratio test
## 
## Model 1: PIRI ~ Treat_Type + avgLD_l + l.Veg_Total + offset(l.TFT) + (1 | 
##     Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + avgLD_l + offset(l.TFT) + (1 | Site/Plot_No)
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1  10 -132.18                     
## 2   9 -132.19 -1 0.0353      0.851
# this works, but ls.m19 has lower AIC
ls.m20 <- glmmTMB(PIRI ~ Treat_Type + avgLD_l + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 ziformula = ~Region,
                 family = nbinom1)
AIC(ls.m20) #286.6
## [1] 286.6293
rm(ls.m13, ls.m14, ls.m15, ls.m16, ls.m17, ls.m18)

Ok, test model fit with treat type (same as model 11)

ls.m19_sr <- simulateResiduals(ls.m19, n = 1000, plot = TRUE) #passes

testResiduals(ls.m19_sr) #pases

## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.037938, p-value = 0.4705
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.0031118, p-value = 0.814
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.00000000 0.01711847
## sample estimates:
## outlier frequency (expected: 0.00112449799196787 ) 
##                                                  0
## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.037938, p-value = 0.4705
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.0031118, p-value = 0.814
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.00000000 0.01711847
## sample estimates:
## outlier frequency (expected: 0.00112449799196787 ) 
##                                                  0
testZeroInflation(ls.m19_sr) #passes

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 0.99996, p-value = 0.972
## alternative hypothesis: two.sided
testQuantiles(ls.m19_sr) #passes

## 
##  Test for location of quantiles via qgam
## 
## data:  simulationOutput
## p-value = 0.3178
## alternative hypothesis: both
#This is a better fit (using nbinom instead of poisson). Convergence issues for nbinom2. Would different distribution change the variables included? - No, I went back and ran model elimination with nbinom1 distro instead of poisson. same variables are significant

ls.m20_sr <- simulateResiduals(ls.m20, n = 1000, plot = TRUE) #passes

testResiduals(ls.m20_sr) #passes

## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.028424, p-value = 0.8159
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.13922, p-value = 0.774
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 1, observations = 498, p-value = 0.48
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.004016064
## sample estimates:
## outlier frequency (expected: 0.00070281124497992 ) 
##                                        0.002008032
## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.028424, p-value = 0.8159
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.13922, p-value = 0.774
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 1, observations = 498, p-value = 0.48
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.004016064
## sample estimates:
## outlier frequency (expected: 0.00070281124497992 ) 
##                                        0.002008032
testZeroInflation(ls.m20_sr) #passes

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 0.98999, p-value = 0.608
## alternative hypothesis: two.sided
testQuantiles(ls.m20_sr) #passes

## 
##  Test for location of quantiles via qgam
## 
## data:  simulationOutput
## p-value = 0.4708
## alternative hypothesis: both

Trying pairwise comparison of model fit:

emmeans(ls.m19, specs = pairwise ~ Treat_Type, adjust = 'Tukey', type = 'response')
## $emmeans
##  Treat_Type response       SE  df asymp.LCL asymp.UCL
##  Control    0.000518 0.000953 Inf  1.41e-05    0.0191
##  FallRx     0.009362 0.016220 Inf  3.14e-04    0.2793
##  Harvest    0.172374 0.280718 Inf  7.08e-03    4.1946
##  MowRx      0.005164 0.008595 Inf  1.98e-04    0.1348
##  SpringRx   0.010190 0.017281 Inf  3.67e-04    0.2829
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast             ratio       SE  df null z.ratio p.value
##  Control / FallRx    0.0553  0.10093 Inf    1  -1.587  0.5060
##  Control / Harvest   0.0030  0.00593 Inf    1  -2.941  0.0271
##  Control / MowRx     0.1003  0.17716 Inf    1  -1.302  0.6901
##  Control / SpringRx  0.0508  0.10026 Inf    1  -1.510  0.5557
##  FallRx / Harvest    0.0543  0.10307 Inf    1  -1.535  0.5395
##  FallRx / MowRx      1.8129  3.06288 Inf    1   0.352  0.9967
##  FallRx / SpringRx   0.9187  1.74321 Inf    1  -0.045  1.0000
##  Harvest / MowRx    33.3770 61.40768 Inf    1   1.907  0.3137
##  Harvest / SpringRx 16.9152 33.56315 Inf    1   1.425  0.6112
##  MowRx / SpringRx    0.5068  0.93049 Inf    1  -0.370  0.9960
## 
## P value adjustment: tukey method for comparing a family of 5 estimates 
## Tests are performed on the log scale

Control vs Harvest is the only significant different (p = 0.0271)

Would things fit better with zi ~Region? Answer: NO, good model fit but higher AIC

GG Plot for ls piri

# PIRI LS -------------------
ls.p1 <- ggpredict(ls.m19, terms = c("avgLD_l", "Treat_Type"))


ggplot(ls.p1, aes(x, predicted, color = group))+
  geom_line(linewidth = 1.25, show.legend = FALSE) +
  scale_color_manual(values = c("#D8B70A",  "#02401B", "#A2A475", "#81A88D", "#972D15"))+
  theme_classic()+
  theme(panel.background = element_blank()) +
  theme(axis.text = element_text(size = 14), axis.title = element_text(size = 20))+
  labs(x = "Average leaf litter depth (log transformed)",
       y = "Pitch pine stems \n (corrected for time from treatment)")

ggsave(filename = "LS_PIRI_avgld.tiff", path = "Plots/PIRI_GLMM", width = 7, height = 5, device = "tiff", dpi = 700)

Shrub oak models

A losing proposition. No models pass fit tests.

Abandon all hope.

Are SO LS present in all regions and treatment types?

tapply(ls.all3$Shrub_Oak, ls.all3$Region, summary)
## $ALB
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    0.00   10.50   18.29   34.00   68.00 
## 
## $LI
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   1.000   4.452   6.000  30.000 
## 
## $MA
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    0.00    8.00   14.49   19.00   73.00 
## 
## $ME
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    8.00   13.00   16.42   23.50   47.00 
## 
## $NH
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    6.25   17.00   19.14   27.00   74.00
tapply(ls.all3$Shrub_Oak, ls.all3$Treat_Type, summary)
## $Control
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   2.000   5.684   9.000  47.000 
## 
## $FallRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    2.25   13.00   16.08   25.00   67.00 
## 
## $Harvest
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   1.500   5.768  10.000  32.000 
## 
## $MowRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00   13.00   26.00   27.68   40.00   74.00 
## 
## $SpringRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     0.0     0.0     7.0    12.6    17.0    61.0
so.ls1 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.BA_piri + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = nbinom1)
AIC(so.ls1) # 3250.6

# test piri ba
so.ls2 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = nbinom1)
AIC(so.ls2) #3251

lrtest(so.ls1, so.ls2) #p = 0.12

# test total BA
so.ls3 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = nbinom1)
AIC(so.ls3) #3251.6

lrtest(so.ls2, so.ls3) # p = 0.1

# test mineral
so.ls4 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = nbinom1)
AIC(so.ls4) #3253.8

lrtest(so.ls3, so.ls4) # p = 0.041

# test avg ld
so.ls5 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = nbinom1)
AIC(so.ls5) #3249.6

lrtest(so.ls5, so.ls3) # p = 0.9

# test piri
so.ls6 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = nbinom1)
AIC(so.ls6) #3250.6

lrtest(so.ls6, so.ls5) # p = 0.08

# test other
so.ls7 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = nbinom1)
AIC(so.ls7) #3256.5

lrtest(so.ls6, so.ls7) # p = 0.005

# test veg
so.ls8 <- glmmTMB(Shrub_Oak ~ Treat_Type +l.other  + l.Mineral + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = nbinom1)
AIC(so.ls8) # 3257.3




lrtest(so.ls8, so.ls6) # p > 0.003
#seems to be best model
so.ls6 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 ziformula = ~1,
                 family = nbinom1)

so.ls6_sr <- simulateResiduals(so.ls6, n=1000, plot = TRUE)

testResiduals(so.ls6_sr)

## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.031295, p-value = 0.7139
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.11972, p-value = 0.002
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.007078313
## sample estimates:
## outlier frequency (expected: 0.000863453815261044 ) 
##                                                   0
## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.031295, p-value = 0.7139
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.11972, p-value = 0.002
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.007078313
## sample estimates:
## outlier frequency (expected: 0.000863453815261044 ) 
##                                                   0
#testing for ziformula: 1, Region + Treat_Type, Treat_Type, Region, l.TFT, combinations therein and with nbinom2 - all are underdispersed. 

Seemingly best models

so.ls5 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = poisson)
AIC(so.ls5) #3378.3 ---- keep mineral, lower AIC

so.ls5_sr <- simulateResiduals(so.ls5, n = 1000, plot = TRUE) #quantile test not looking good
testResiduals(so.ls5_sr)
testDispersion(so.ls5_sr, alternative = "less") #again, under dispersed

so.ls5a <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = compois)
AIC(so.ls5a) #3290

so.ls5a_sr <- simulateResiduals(so.ls5a, n = 1000, plot = TRUE) 
testResiduals(so.ls5a_sr) #still underdispersed

so.ls5b <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = genpois)
AIC(so.ls5b) #3274.5

so.ls5b_sr <- simulateResiduals(so.ls5b, n = 1000, plot = TRUE)
testResiduals(so.ls5b_sr) #still underdispersed

Test second model

so.ls6 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + l.BA_HA + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = poisson)
AIC(so.ls6) #3380

so.ls6_sr <- simulateResiduals(so.ls6, n = 1000, plot = TRUE)
testResiduals(so.ls6_sr) # fails dispersion

so.ls6a <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + l.BA_HA + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = genpois)
AIC(so.ls6a) #3275.5
#compois distro fails and won't produce DHARMa results
so.ls6a_sr <- simulateResiduals(so.ls6a, n = 1000, plot = T)
testResiduals(so.ls6a_sr) # fails dispersion

testing other distros

so.ls6 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + l.BA_HA + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = poisson)
AIC(so.ls6)

so.ls6a <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + l.BA_HA + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = genpois)
AIC(so.ls6a)

# compois distro fails, won't produce DHARMa results

so.ls6b <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + l.BA_HA + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = tweedie)
AIC(so.ls6b)

so.ls6b_sr <- simulateResiduals(so.ls6b, n = 1000, plot = TRUE)
testDispersion(so.ls6b_sr) #fails

so.ls6c <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + l.BA_HA + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = t_family)
AIC(so.ls6c)

so.ls6c_sr <- simulateResiduals(so.ls6c, n = 1000, plot = TRUE)
testResiduals(so.ls6c_sr)
testDispersion(so.ls6c_sr) #passes
testQuantiles(so.ls6c_sr)


so.ls6d <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + l.BA_HA + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = gaussian)
AIC(so.ls6d)

so.ls6d_sr <- simulateResiduals(so.ls6d, n = 1000, plot = TRUE) #fails
plot(so.ls6d_sr)
testResiduals(so.ls6d_sr)

Pairwise with treatment type

emmeans(so.ls6, specs = pairwise ~ Treat_Type, adjust = 'Tukey', type = 'response')
hist(ls.all3$Shrub_Oak)

ggplot(ls.all3, aes(x=Shrub_Oak))+
  geom_histogram(binwidth = 2)+
  facet_grid(rows = vars(Treat_Type))

Graph shrub oak models

# LS SO plot 1 ------------------- other vs. TT
plot_model(so.ls6, 
                       type = 'pred', 
                       terms = c('l.other', 'Treat_Type'),
                       axis.title = c("Other small seedling counts (log transformed)", "Total Count of Shrub Oak"),
                       title = "Predicted Counts of Shrub Oak \n Seedlings >/=50cm & <2.5cm DBH",
           legend.title = "Treatment Type",
           line.size = 1,
           value.offset = 'Treat_Type',
           ci.lvl = NA,
           colors = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15")) 

# LS SO plot 2 ------------------- BA_10m vs. TT
plot_model(so.ls6, 
                       type = 'pred', 
                       terms = c('l.other', 'Treat_Type'),
                       axis.title = c("Total Basal Area per 10 square meters (log transformed)", "Total Count of Shrub Oak"),
                       title = "Predicted Counts of Shrub Oak \n Seedlings >/=50cm & <2.5cm DBH",
           legend.title = "Treatment Type",
           line.size = 1,
           value.offset = 'Treat_Type',
           ci.lvl = NA,
           colors = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15")) 

# LS SO plot 3 ------------------- l.mineral vs. TT
plot_model(so.ls6, 
                       type = 'pred', 
                       terms = c('l.other', 'Treat_Type'),
                       axis.title = c("Amount of exposed mineral soil (log transformed)", "Total Count of Shrub Oak"),
                       title = "Predicted Counts of Shrub Oak \n Seedlings >/=50cm & <2.5cm DBH",
           legend.title = "Treatment Type",
           line.size = 1,
           value.offset = 'Treat_Type',
           ci.lvl = NA,
           colors = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15")) 

#Start model fit with genpois distro and then do elimination?

so.ls1 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.BA_piri + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = genpois)
AIC(so.ls1) #3275.5


so.ls1b <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.BA_piri + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = compois)
AIC(so.ls1b) #this is taking a very long time - so i quit, moving forward with genpois

#test piri ba
so.ls2 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = genpois)
AIC(so.ls2) #3276.4

lrtest(so.ls1, so.ls2) #drop

#test ba
so.ls3 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = genpois)
AIC(so.ls3) #3276.8

#test mineral
so.ls4 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = genpois)
AIC(so.ls4) #3279.2

lrtest(so.ls3, so.ls4) #keep mineral

# test ld
so.ls5 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = genpois)
AIC(so.ls5) #3274.8

# test other
so.ls6 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI  + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = genpois)
AIC(so.ls6) #3280

lrtest(so.ls6, so.ls5) #keep other

# test piri
so.ls7 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = genpois)
AIC(so.ls7) #3275.9

lrtest(so.ls7, so.ls5) #drop piri

#going to drop veg due to collinearity issues
# test ld
so.ls8 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + l.Mineral + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = genpois)
AIC(so.ls8) #3284

so.ls8_sr <- simulateResiduals(so.ls8, n = 1000, plot = TRUE)
testResiduals(so.ls8_sr)
testDispersion(so.ls8_sr, alternative = "less")
#underdispersed
testZeroInflation(so.ls8_sr)

so.ls9 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + l.Mineral + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = compois)
AIC(so.ls9) #3301.9

so.ls9_sr <- simulateResiduals(so.ls9, n = 1000, plot = TRUE)
testResiduals(so.ls9_sr) #worse on dispersion

so.ls10 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + l.Mineral + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = nbinom1)
AIC(so.ls10) #3257.3

so.ls10_sr <- simulateResiduals(so.ls10, n = 1000, plot = TRUE)
testResiduals(so.ls10_sr) #fails dispersion


so.ls11 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + l.Mineral + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = nbinom2)
AIC(so.ls11) #3352

so.ls11_sr <- simulateResiduals(so.ls11, n = 1000, plot = TRUE)
testResiduals(so.ls11_sr) # fails dispersion
testZeroInflation(so.ls11_sr)

# run poisson model with zero inflation?
so.ls12 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + l.Mineral + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 ziformula = ~1,
                 family = genpois)
AIC(so.ls12) #3271.5

so.ls12_sr <- simulateResiduals(so.ls12, n = 1000, plot = TRUE)
testResiduals(so.ls12_sr) #still fails dispersion
testZeroInflation(so.ls12_sr)

so.ls13 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + l.Mineral + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = poisson)
AIC(so.ls13)

so.ls13_sr <- simulateResiduals(so.ls13, n = 1000, plot = TRUE)
testResiduals(so.ls13_sr) #still fails dispersion
testZeroInflation(so.ls13_sr) #does fail

ok, zero inflation issues

#starting with veg dropped due to collinearity issues
so.ls1 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.BA_piri + l.Mineral + avgLD_l + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 ziformula = ~1,
                 family = nbinom1)
AIC(so.ls1) #3253.5

#test piri ba
so.ls2 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.Mineral + avgLD_l + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 ziformula = ~1,
                 family = nbinom1)
AIC(so.ls2) #3252

#test ba
so.ls3 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.Mineral + avgLD_l + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 ziformula = ~1,
                 family = nbinom1)
AIC(so.ls3) #3254

lrtest(so.ls2, so.ls3) # p = 0.6, drop

#test ld
so.ls4 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.Mineral + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 ziformula = ~1,
                 family = nbinom1)
AIC(so.ls4) #3252

#test piri
so.ls5 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + l.Mineral + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 ziformula = ~1,
                 family = nbinom1)
AIC(so.ls5) #3254

lrtest(so.ls4, so.ls5) # keep piri

#test other
so.ls6 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI  + l.Mineral + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 ziformula = ~1,
                 family = nbinom1)
AIC(so.ls6) # 3256

lrtest(so.ls4, so.ls6) #keep other

#test mineral
so.ls7 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 ziformula = ~1,
                 family = nbinom1)
AIC(so.ls7) #3255.4

lrtest(so.ls7, so.ls4) #keep mineral

#test treat type
so.ls8 <- glmmTMB(Shrub_Oak ~ l.PIRI + l.other + l.Mineral + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 ziformula = ~1,
                 family = nbinom1)
AIC(so.ls8) #3341.4

lrtest(so.ls4, so.ls8) # keep tt

seemingly best model

so.ls4 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.Mineral + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 ziformula = ~1,
                 family = nbinom1)

so.ls4_sr <- simulateResiduals(so.ls4, n = 1000, plot = TRUE)
testZeroInflation(so.ls4_sr) #passes
testResiduals(so.ls4_sr) #fails dispersion

so.ls4b <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.Mineral + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 ziformula = ~Treat_Type,
                 family = nbinom1)

so.ls4b_sr <- simulateResiduals(so.ls4b, n = 1000, plot = TRUE)
testZeroInflation(so.ls4b_sr) #passes
testResiduals(so.ls4b_sr)

#try with second nbinom distro

so.ls4c <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.Mineral + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 ziformula = ~1,
                 family = nbinom2)

so.ls4c_sr <- simulateResiduals(so.ls4c, n = 1000, plot = TRUE)
testZeroInflation(so.ls4c_sr) #passes
testResiduals(so.ls4c_sr) #worse

so.ls4d <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.Mineral + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 ziformula = ~Treat_Type,
                 family = nbinom2)

so.ls4d_sr <- simulateResiduals(so.ls4d, n = 1000, plot = TRUE)
testZeroInflation(so.ls4d_sr) #passes
testResiduals(so.ls4d_sr) #worse

# maybe 
so.ls4e <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.Mineral + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 ziformula = ~Treat_Type,
                 family = genpois)

so.ls4e_sr <- simulateResiduals(so.ls4e, n = 1000, plot = TRUE)
testZeroInflation(so.ls4e_sr) #passes
testResiduals(so.ls4e_sr) #fails again

starting modeling again. going to model with poisson and check out zero inflation issues. then reduce model or try other distros

m1 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.Mineral + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = poisson)
AIC(m1) #3391.3

m1_sr <- simulateResiduals(m1, n = 1000, plot = T)
testZeroInflation(m1_sr)#fails

m2 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.Mineral + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
              ziformula = ~1,
                 family = poisson)
#1, region, and tt all fail with poisson distro

m2 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.Mineral + offset(l.TFT) + (1|Site/Plot_No),
              data = ls.all3,
              family = nbinom1)
AIC(m2) #3254.1

m2_sr <- simulateResiduals(m2, n = 1000, plot = T)
testZeroInflation(m2_sr) #passes
testResiduals(m2_sr) # fails, is underdispersed
testDispersion(m2_sr, alternative = "less") # = 0.006

#try nbimom2
m3 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.Mineral + offset(l.TFT) + (1|Site/Plot_No),
              data = ls.all3,
              family = nbinom2)
AIC(m3) #3313.1

m3_sr <- simulateResiduals(m3, n = 1000, plot = T)
testZeroInflation(m3_sr) #fails
testResiduals(m3_sr) # fails, is underdispersed
testDispersion(m3_sr, alternative = "less") # <0.001

#try genpois (it fails)
#try compois - taking a very long time to run
m4 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.Mineral + offset(l.TFT) + (1|Site/Plot_No),
              data = ls.all3,
              ziformula = ~1,
              family = compois)
AIC(m4) #3270

m4_sr <- simulateResiduals(m4, n = 1000, plot = T)
testZeroInflation(m4_sr) #passes
testResiduals(m4_sr) # fails, is underdispersed
testDispersion(m4_sr, alternative = "less") # <0.001

#compois is worse

m5 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.Mineral + offset(l.TFT) + (1|Site/Plot_No),
              data = ls.all3,
              family = nbinom1)

m5_sr <- simulateResiduals(m5, n = 1000, plot = T)
testZeroInflation(m5_sr) #passes
testResiduals(m5_sr) #0.006

#region is slightly worse (0.002) for ziformula; treat_type is 0.006; site fails, l.TFT is 0.006;BAHA fails 

# maybe having less variables will change the sr model

m6 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.Mineral + offset(l.TFT) + (1|Site/Plot_No),
              data = ls.all3,
              family = nbinom1)

m6_sr <- simulateResiduals(m6, n = 1000, plot = T)
testZeroInflation(m6_sr) #passes
testResiduals(m6_sr) #worse 0.002

m7 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + l.Mineral + offset(l.TFT) + (1|Site/Plot_No),
              data = ls.all3,
              ziformula = ~1,
              family = nbinom1)

m7_sr <- simulateResiduals(m7, n = 1000, plot = T)
testZeroInflation(m7_sr) #passes
testResiduals(m7_sr) #0.008


m8 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.Mineral + offset(l.TFT) + (1|Site/Plot_No),
              data = ls.all3,
              ziformula = ~1,
              family = nbinom1)

m8_sr <- simulateResiduals(m8, n = 1000, plot = T)
testZeroInflation(m8_sr) #passes
testResiduals(m8_sr) #0.002

m10 <- glmmTMB(Shrub_Oak ~ Treat_Type + offset(l.TFT) + (1|Site/Plot_No),
              data = ls.all3,
              ziformula = ~1,
              family = nbinom1)

m10_sr <- simulateResiduals(m10, n = 1000, plot = T)
testZeroInflation(m10_sr) #passes
testResiduals(m10_sr) #0.004

m11 <- glmmTMB(Shrub_Oak ~ Treat_Type + offset(l.TFT) + (1|Site/Plot_No),
              data = ls.all3,
              family = nbinom1)

m11_sr <- simulateResiduals(m11, n = 1000, plot = T)
testZeroInflation(m11_sr) #passes
testResiduals(m11_sr) #0.004


m11 <- glmmTMB(Shrub_Oak ~ Treat_Type + offset(l.TFT) + (1|Site/Plot_No),
              data = ls.all3,
              ziformula = ~1,
              family = gaussian)

m11_sr <- simulateResiduals(m11, n = 1000, plot = T)
testZeroInflation(m11_sr) #passes
testResiduals(m11_sr) #0.004


m2 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.Mineral + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
              ziformula = ~1,
                 family = gaussian)
AIC(m2)
m2_sr <- simulateResiduals(m2, n = 1000, plot = T)
testZeroInflation(m2)
testResiduals(m2) #fails uniformity, passes dispersion

m2 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.Mineral + offset(l.TFT) + (1|Site/Plot_No),
              data = ls.all3,
              ziformula = ~1,
              family = tweedie)
AIC(m2)
m2_sr <- simulateResiduals(m2, n = 1000, plot = T)
testZeroInflation(m2)
testResiduals(m2) #fails dispersion passes uniformity

m2 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.Mineral + offset(l.TFT) + (1|Site/Plot_No),
              data = ls.all3,
              family = tweedie)
AIC(m2)
m2_sr <- simulateResiduals(m2, n = 1000, plot = T)
testZeroInflation(m2)
testResiduals(m2) #fails dispersion passes uniformity

nbinom distro 1 solves 0 inflation issues, but models are still underdispersed

I’ve tried poisson, nbinom1, nbinom1 with ZI, nbinom2, genpois, compois, tweedie, gaussian and none of these distros work. They are all underdispersed

These are very full models with many different distributions. All are underdispersed.

m1 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.BA_piri + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = nbinom1)
AIC(m1) 
## [1] 3250.567
m1_sr <- simulateResiduals(m1, n = 1000, plot = T)

testZeroInflation(m1_sr) #passes

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0085, p-value = 0.96
## alternative hypothesis: two.sided
testResiduals(m1_sr) #under dispersed, 0.002

## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.030372, p-value = 0.7479
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.18917, p-value = 0.006
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.006024096
## sample estimates:
## outlier frequency (expected: 0.00106425702811245 ) 
##                                                  0
## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.030372, p-value = 0.7479
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.18917, p-value = 0.006
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.006024096
## sample estimates:
## outlier frequency (expected: 0.00106425702811245 ) 
##                                                  0
m2 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.BA_piri + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = genpois)
AIC(m2) 
## [1] 3275.482
m2_sr <- simulateResiduals(m2, n = 1000, plot = T)

testZeroInflation(m2_sr) #passes

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.1026, p-value = 0.49
## alternative hypothesis: two.sided
testResiduals(m2_sr) #under dispersed, 0.004

## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.040086, p-value = 0.4003
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.12405, p-value = 0.004
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.008032129
## sample estimates:
## outlier frequency (expected: 0.00162650602409639 ) 
##                                                  0
## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.040086, p-value = 0.4003
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.12405, p-value = 0.004
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.008032129
## sample estimates:
## outlier frequency (expected: 0.00162650602409639 ) 
##                                                  0
m3 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.BA_piri + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = ls.all3,
                 family = compois)
AIC(m3) #this is taking a very long time to run
## [1] 3291.162
m3_sr <- simulateResiduals(m3, n = 1000, plot = T)

testZeroInflation(m3_sr) #fails

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.3471, p-value = 0.038
## alternative hypothesis: two.sided
testResiduals(m3_sr) #fails dispersion

## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.025493, p-value = 0.9026
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.045771, p-value = 0.002
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.00000000 0.01310241
## sample estimates:
## outlier frequency (expected: 0.00152610441767068 ) 
##                                                  0
## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.025493, p-value = 0.9026
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.045771, p-value = 0.002
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.00000000 0.01310241
## sample estimates:
## outlier frequency (expected: 0.00152610441767068 ) 
##                                                  0
m4 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.Mineral +  offset(l.TFT) + (1|Site/Plot_No),
              data = ls.all3,
              ziformula = ~1,
              family = compois)
AIC(m4)
## [1] 3269.97
m4_sr <- simulateResiduals(m4, n = 1000, plot = T)

testZeroInflation(m4_sr) #passes

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.097, p-value = 0.462
## alternative hypothesis: two.sided
testResiduals(m4_sr) #fails dispersion

## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.028524, p-value = 0.8125
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.051874, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.00000000 0.01817269
## sample estimates:
## outlier frequency (expected: 0.00144578313253012 ) 
##                                                  0
## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.028524, p-value = 0.8125
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.051874, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.00000000 0.01817269
## sample estimates:
## outlier frequency (expected: 0.00144578313253012 ) 
##                                                  0

I cannot find a model that works.